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ABSTRACT: 

The equations governing free turbulent mixing are derived from the 
Navier -Stokes equations and transformed into a mathematical plane which 
is explicitly independent of the eddy viscosity model. The coupled 
momentum and turbulent kinetic energy equations are analytically solved 
in the transformed plane by a perturbation technique and subsequently 
retransformed into physical space based on a hypothesized dependence of 
the eddy viscosity on the turbulent kinetic energy. The adequacy of a 
given model in reproducing the velocity and turbulent kinetic energy 
field is assessed by comparing the results of the analysis with some 
experimental data of planar turbulent wake mixing in constant adverse 
and favorable pressure gradients. 
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NOMENCLATURE 



Cross-Sectional Area of Wind Tunnel Test Section 
Strength of Pressure Gradient, Equations (49) and (51) 
Constant in Equation (32) 

Turbulent Dissipation 
Displacement Thickness 

Perturbation Parameter = 1/u q ; DK/p in Equation 24 

Arbitrary Function in Equation (73) 

Fe -Jf( c p)dcp 

Greens Function 

Shape Factor = 6*/9 

Initial Conditions in Equation (75) 

Unit Vector in x-direction 
Nonhomogene ous Function in 
Unit Vector in y-direction 
Turbulent Kinetic Energy = 

Constant in Equation (52) 

Unit Vector in z-direction 

2 -2 
u - u 
e 

Macroscale of Turbulence 
Coefficient of Viscosity 
Eddy Viscosity 
Kinematic Viscosity 
Turbulent Kinematic Viscosity 
Pressure 

Stream Function Defined in Equation (25) 

Transformed Independent Variable in Equation (33) 
Dynamic Pressure 



Equation (74) = g(cp>Y)e~/ f 



/ 2 2 2 \ 
\ \u' + v' + w' ) 



iv 



p 


Density 


°K = 


Turbulent Diffusion Coefficient in Equation (23) 


S 


2 

Velocity Ratio Function = (u q /u ) 


t 


Time 


T 


Apparent Turbulent Shear Stress 


e 


Momentum Thickness 


u = 


Axial Velocity (x-direction) 


U 

e 


Average Axial Velocity External to the Wake 


V 


Lateral Velocity (y-direction) 


V 


Velocity Vector = ui + vj + wk 


w 


Transverse Velocity (z-direction) 


X 


Axial Coordinate 


y 


Lateral Coordinate 


z = 


Transverse Coordinate 



Subscripts 



e = 


External to the Wake 


= 


Three Orthogonal Components of a Vector 


L 


Exit of Test Section 


0 


Inlet of Test Section 


t 


Turbulent Function 


0,1 = 


Zeroth and First Order Solutions of x and K 



Superscripts 



u = 


Time Average of u 


u' 


Fluctuating Part of u 



v 



INTRODUCTION 



The solution of the Navier-Stokes equations as applied to the problem 
of turbulent flow has classically been approached either from the point of 
view of elegant rigor on limit ingly simple flows or one of empirically 
guided analysis of flows of real interest. The Reynolds averaging 
technique is nearly a pre-requisite to attacking any real turbulent flow 
and the associated necessity to empirically close the system of equations 
with a Reynolds stress model makes the mathematics tractable yet often 
unreliable. Various so-called eddy viscosities which arise from the 
Boussinesq laminarizat ion of the apparent turbulent stresses are functions 
of the flow field and as yet there are no known universal functions which 
adequately model these stresses in all cases. 

Possibly the most successful application of eddy viscosity approaches 
lies in the area of free mixing where velocity differences through the flow 
field are small and the associated turbulent field is a fairly simple one. 
However, in the pressure gradient situations which are encountered in 
ejector and combustor mixing phenomena, classical eddy viscosities fail 
since they are explicitly independent of the turbulent field which, in 
these cases, can become complex and exert a dominant effect on the form of 
the eddy viscosity. Modern approaches have attempted to include the local 
turbulence structure effect on the eddy viscosity through an explicit 
dependence on the local turbulent kinetic energy. 

In conjunction with a highly idealized wake mixing experiment, the 
two-dimensional incompressible turbulent wakes in constant adverse and 
favorable pressure gradients have been studied analytically with a formu- 
lation of the equations which allows the coordinated solutions of the 
velocity and turbulence field in a transformed plane which is explicitly 
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independent of the eddy viscosity model. Subsequently, the retransform- 
ation of the equations for a variety of models allows for a comparison 
of the adequacy of the models and an evaluation as to which most 
accurately reproduces both the velocity and kinetic energy fields. 



The basis for any rigorous analysis of problems involving turbulent 
fluid flow must, to the best of current knowledge, be founded in the 
Navier-Stokes equations. For an incompressible, Newtonian fluid these 
may be written 



No analytical solution of the full equations appears possible and, 
although some success has been shown with numerical approaches to the 
solution of the equations, as of yet no numerical attack for a genuinely 
turbulent flow is feasible (Ref. l). Apart from the numerical approach, 
only Fourier analysis has shown any progress in the solution of the 
equations. However, with this method, only limit ingly simple turbulent 
fields have been treated and its relevance to a general turbulent mixing 
problem has yet to be proven. 

Classically, the most fruitful approach in the analysis of turbulence 
has been to decompose each of the dependent variables in Equations 1 and 
2 into a mean term, which is independent of time or has a long character- 
istic period with respect to the turbulent fluctuations, plus a fluctua- 
ting term whose time average is zero. The velocity field 



EQUATION DEVELOPMENT 



V • V = 0 



(i) 




( 2 ) 



V = ui + vj + wk 



( 3 ) 
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may be decomposed term by term such as the decomposition of u 



u (x,y,z,t) = u (x,y,z) + u* (x,y,z,t) 



( 4 ) 



When this decomposition is applied to each term in the continuity 
equation we obtain 



du + Sv + dw + du' + dv 1 + dw 1 
dx dy dz dx dy dz 



( 5 ) 



If we now take the time average of Equation ( 5 ), we obtain 



Bu + dv + dw = o 
dx dy dz 



( 6 ) 



which is the continuity equation for the mean flow. Upon subtracting 
Equation (6) from Equation (5), we obtain 



du' 

dx 



+ 





( 7 ) 



This is the continuity equation which the velocity fluctuations must 
satisfy. Prior to applying this decomposition technique to the momentum 
equations we first reformulate the convective operator in a conservative 
form, with the aid of the continuity equation. The general convective 
operator 



_d_ 

dt 



+ u 



_d_ 

dx 



+ v 



d_ 

By 



+ w 



_d_ 

dz 



( 8 ) 



becomes 



_d_ 

dt 



+ 



_d_ 

dx 



u < > >♦&»< > 



( 9 ) 
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where the appropriate dependent variable is placed within the parentheses . 
Applying this formulation to Equations (2) and expanding them in rectang- 
ular Cartesian coordinates, we obtain 



|jr + £- (u 2 ) + (uv) + (uw) = 
St 5x Sy Sz 



Up , 2 

— r*- + v V u 
p Sx 



TT + (uv) + (v 2 ) + (wv) 

St Sx Sy Sz 



- i |£ - v V 2 v 
p Sy 



(10) 



+ A" (uw) + (vw) + (w 2 ) 

St Sx Sy Sz 



1 Sp ^ 2 

= - + \. V w 



P Sz 



Each of the dependent variables in Equations (10) is now decomposed and 
the appropriate expression is inserted into the equations. The resulting 
system of equations is 



Su + Su 1 + Su 2 
St St Sx 



Suu' Su' 2 Suv 

Sx Sx Sy 



+ 



Suv' Su 1 v Su' v' 

Sy Sy Sy 



Suw + Su'w + Suw' 
Sz Sz Sz 



+ 5u'w ’ 
Sz 



1 SP 
p Sx 



1921 + 

p Sx 



„2- A 

v V u + 



V v u' 



Sv + 5 v' + Suv + Su' v + Suv* + Su' v' + Sv 2 + 2 Sw' + Sv' 2 

St st Sx Sx Sx Sx Sy Sy Sy 

( 11 ) 



Swv Sw' v + Swv* + Sw'v' 
Sz Sz Sz Sz 



1 Sp 

p 9y 



1 SP' , 
P Sy 



2 - 

v V v + v 



2 

V v 



I 
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dw + 8w' + puw du’w ouw' du'w 1 + dvw + dvw' + 3v'w dv'w' 

dt dt dx dx dx dx dy dy dy dy 



+ 




_ 2 

dww 1 dw' 
dz dz 



1 dp 
p d z 



1 dpj_ 
p dz 



2 — 2 
+V^w+ V vw' 



If we now take the time average of each equation, we can simplify the 
system by dropping out all those terms whose time average is zero (i.e. 
those which have terms that are linear in the fluctuating properties). 
Note, however, that we must retain the nonlinear terms containing 
fluctuations since the products or powers of purely fluctuating terms may 
generate a steady time averaged value. When the appropriate time averages 
are taken of each term in Equations (ll) we retain the following system of 
equations 



+ g- (u 2 ) + (w) + g- (w) = 

dt dx dy dz 



- — (p + pu' 2 ) - u'v' - u'w' 



P dx 



dy 



dz 



2 - 

+ vi V u 



§ + 1 <-> + & + h <-> - 



2 - 



- j k (p + pv ' 2) ' & U ' T ' • h ' ,w ' + v v_v 



dw 3 / \ d / \ , d / — 2 n 1 d / |2\ 

3t + (uw) + 37 (w) SI ( " J = ' P SI (p p ” ) 



u'w' 

dx 



v'w' + \j V 



dy 



2 - 

w 



( 12 ) 



Note that the price of this "simplification" through the Reynolds averaging 
has been the introduction of six new unknowns by discarding of all of the 
phase information of the fluctuations. For a steady, two-dimensional mean 
flow Equations (12) may be simplified to 



+ Sv = o 
dx dy 



( 13 ) 
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( 14 ) 





2 

In general we may make the assumption that the apparent pressures (pu 1 , 
pv 1 ) may be neglected with respect to p, We should note that when p = 



constant some care must be exercised in applying this approximation since 



for a particular problem. In addition, for a "thin 11 wake, we may make the 
standard boundary layer-type parallel flow approximation that 



This eliminates the necessity of solving the lateral momentum equation, 
since the only unknowns are u and v with p being imposed on the mixing 
region by the external flow. The only remaining term which explicitly 
involves the turbulent field is -u ! v ! which represents an apparent Reynolds 
shearing stress (t). 



In general this apparent stress is very large with respect to the average 
laminar shear stress and we may neglect the laminar component. With these 
approximations the resulting system of equations to be solved is 





(16) 



T = “ pU ! V ! 



( 17 ) 




dx 



ay 



(18) 




1 dp + 1 dr 

p dx p 3y 



( 19 ) 
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At this point, some empiricism is necessary since t is not retrievable 
from the equation set we have derived. This is clear from Equations (ll) 
when either is multiplied by the appropriate fluctuating velocity and 
time averaging is performed to obtain an equation for the Reynolds stress, 
triple correlations appear which axe unknowns also. This cascade of 
unknown higher order correlations continues for all further equations 
which axe derived. Some success (Ref. 2) has been made attacking these 
problems by making purely heuristic approximations in higher order 
correlations in an effort to make the least sensitive approximations 
possible in the equations. However, in lieu of attacking this spiralling 
set of equations, it has generally proved more effective to introduce some 
empirically based models of the Reynolds stress into Equations (l8) and 
(19). The most successful models have been based on extensions of 
Prandtl's mixing length analysis which analogizes turbulent eddy momentum 
transfer with the molecular manifestation of viscosity. In line with this 
approach an eddy viscosity, is introduced into the problem via the 
definition 



where is a function of the local mean flow field. Classical models 



however, more modern results (Ref. 3 ) indicate that some specific 
dependence on the turbulence field is indicated. Thus we hypothesize that 




( 20 ) 



infer that is purely a function of the local non-turbulent mean flow, 



|J. t = (u,v,K) 



where 



K s \ (u' 2 + v' 2 + w' 2 ) 



( 21 ) 
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In order to implement a model of this form, we must formulate and solve an 



equation for K along with the momentum and continuity equations. In order 
to obtain an equation for K, we multiply Equations (ll) by u' , v' , w' 
respectively, add them and take the time average. This operation results 
in the following equation for K which is most conveniently written in 
tensor notation 



DK 

Dt 



_c>_ 

ax. 



f— u'.p' +^-u'u'.u'. - v u ' - ( 
\pj 2 k 1 j 1 \ 



^i 

d x. 



au* . 



Sx 



t)) 



(22) 



u .u 



if 



au . 



Vdx, 



Su 



V 



(t^ + 

u 
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The respective convection, diffusion, production, and dissipation terms 
have been modelled by Patankar and Spalding (Ref. 4) in the following 
equation 



- SK - 3K _ 9 T v t 9K*1 . /9u\ 2 D K 

u 3x v 9y 5y La K SyJ v t \Sy/ ” p 



(23) 



Where a v is the effective Prandtl number for the diffusion of K and D is 

the turbulent dissipation. With appropriate auxiliary expressions or 

equations for a > D , and o we may consider the system of equations to 
t K K 

be closed. Clearly the adequacy of the equation system in modelling any 
particular flow field is dependent on the exact formulation which is used 
for the unknown coefficients . 

Based on dimensional reasoning, we can specify the coefficients 
and D v to be 

IS. 

= const, x density x velocity x length 

o 

D^ = const, x density x velocity” 1 / length 
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From experimental results (Ref. 3) we specify a to be a purely empirical 

ft 



constant in the range 0.5 — 1.0. 

With the explicit expression which we will test for jj, , along with 

T> 

a value for a , only the formulation for D„ is left unspecified. Jones 
ft ft 

(Ref. 5) has hypothesized an equation for e = D^/p of the following form 



De _ 5 / v t d e \ ~ e_ ^"i 

Dt dx. \ct dx .) 1 K v t dx. 

J e j 



du_. du_. 

I — 



. \dx. 
J i 



du 

dx 






0 £ 

K 



(24) 



where and c 2 are empirical constants. In lieu of this complication, 

we will make some experimentally justified approximations for and its 

ft 

dependence on K and the mean velocity field. 

MATHEMATICAL ANALYSIS 

I. Equation Reformulation 

In order to reduce Equations (l8) , (19)? and (23) to a form more 
amenable to an analytic approach, we introduce the stream function 

Y = u , y = - v (25) 

J X* 

which automatically satisfies the continuity equation. With this definition 
of the stream function, the von Mises transformation is applied to the 
independent variables (x,y) to transform them into (x,y) via the following 
equations 



d_ = _d_ 
dx dx 




_d_ = - d_ 
dy U dY 



(26) 

(26) 



9 



When this transformation is applied to Equations (19) and (23), we obtain 



du _ _ _1_ dp + _d_ / _ Bu\ 

dx pu dx v' t U dy/ 



(27) 



dK 

dx 




+ 





(28) 



where continuity need not be solved since v does not appear explicitly in 
the equations when u and K are expressed in stream function variables. 

The dependent variable u is now transformed into X via the equation 

X = u e 2 - u 2 (29) 

2 i 

Upon substituting u = (u - x) 2 into Equations (27) and (28) we can write 
the equations of momentum and turbulent kinetic energy as 

I - «. (1 - fk (30) 

u e ^ 



dK 

dx 



v t U e 







K 



H.u 



t e 




2 } (3D 



To this point, based on the assumptions previously outlined, these 
equations are as exact as the specifications of D^, v^, and To 

consistent with our previous dimensional reasoning, the coefficient 
the kinetic energy equation may be expressed dimensionally as 



be 

H 



in 



— = constant x (velocity/length) 
^t 
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Following experimental indications we hypothesize the exact formulation 
for this coefficient to be 



H t = C x K/£ 2 (32) 

where i is the maxroscale of the turbulent flow, generally accepted to be 
of the order of the mixing width and a physical measure of those turbulent 
eddies most intimately involved in momentum transfer. We have not yet 
explicitly used an independent expression for in the equations and, as 
we shall see, the hypothesis for its formulation need not be specified 
until the equations have been solved in a transformed plane. 

With these equations and expressions, we again transform the indepen- 
dent variables, this time from (x,y) to (cp , Y ) using the following defini- 
tion of cp 

x 

cp = f u v+ dx (33) 

e t 

where we have introduced a crucial yet common and experimentally justified 

assumption that v. = v + ( x ) alone. However, our interest lies in specifying 

a useful v^.(k) which is implicitly [ K (cp )] in the transformed plane and 

we are free to select a particular f . for optimum results, thus v + (x) is in 

J t 

actuality = v t [K(cp,Yj)]. 

Applying this transformation to the equations, we can write the 
momentum and turbulent kinetic energy equations as 



SX 

Sep 



1 - 



I .2 

JO d ■ 



u 



\ <LA 

' By 2 



(34) 
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aK 

Bcp 



i 



III 


17 1 _ JL_\ 2 dK~ 


. i / 


ct k dY 


1 ^ 
ISO 

c\j" 

1 

2 

i 

H 


C\J 

i 



e e 




2 2 
£ U e 




(35) 



Subject to the following boundary and initial conditions: 
B.C.x-*0;K-*0 as y - = “ 

I.C. X = X i (Y) ; K = K.(y) at cp = o 



(36) 



II. Equation Simplification — Linearized Case 

For many wake or jet -like parallel flows the velocity within the wake 
differs only slightly from that in the local external flow. In these cases 



X 



2 

u 

e 



« 1 



and we can write Equations (3*0 and (35) as 



62 L = 9JL 
9cp Sy 2 



(37) 



dK _1_ dfk _ / C \ K 1 /dX\ 2 

^ CT K 9y 2 Vu^ 4u e 2 



(38) 



subject to the identical boundary and initial conditions specified for 
Equations (34) and (35). 



III. Equation Simplif ication--Nonlinearized Case 

In free mixing cases, although the square of the velocity defect (^) 

2 

is small with respect to u , it is not negligible and some influence of 
the finiteness of the velocity defect must be included. With the governing 
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equations in the (cp,Y) plane it is particularly convenient to include 
this effect. This is accomplished, first with the momentum equation, by 
expanding the dependent variables in powers of the parameter which makes 
X small, namely some measure of the freestream velocity. 

With the following definitions 



S 




1 




the momentum equation becomes 




( 39 ) 



( 40 ) 



We are seeking here only a first order correction to the linearized set of 
equations, however, successive higher order approximations may be derived 
in exactly the same manner. Using the binomial theorem, we expand the 
diffusive term and, keeping only a first order correction to the linearized 
case, we can write the momentum equation 

& - (1 - i . S X ) W 

* By 

00 

We now expand the dependent variable x i n a power series x = T, e X- 

i=0 1 

and, in line with the objective of obtaining a first order correction to 
the linearized system, we retain only the first two terms in the series 
which results in 



X = X 0 + e X x 
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Upon substituting this expression into Equation (4l) and collecting terms 
of equal order in e we obtain the following equations for X-^ 



9Xq 

dcp 




(42) 



3Xi 

Sep 



M 2 



i s 






s 2 x 0 

, 2 
dY 



(43) 



Clearly the zeroth order equation is the strictly linearized case and the 
first order solution is the correction term for the finiteness of the 
velocity defect which tends to force the solution to satisfy the full 
equation. 

Similarly, the linearized solution for K may he corrected for the 
improved velocity field solution and for the direct effect of a finite 
velocity defect. Here again we seek a first order correction for the 
linearized solution by expanding ^ and K in Equation (35) in terms of e. 

The expressions for u are expanded with the binomial theorem with use of 
the appropriate expressions for x = X 0 + e the appropriate zeroth and 
first order solutions to the momentum equation. Clearly the equation for 
K is now explicitly dependent upon e and we can expand K also in powers of 
e , again to first order, so that a correction to the linearized case may be 
obtained. Upon substituting 




+ e K 



1 



the respective equations for K Q and can be obtained by collecting terms 
of equal order in e . The equations for K Q and K can be written 
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*0 ±*X + i A ox 

a K 3 t 2 'i 2 u 2 0 ku 2 ^' l) 

e e 



(44) 



^1 _ _1_ 5 K i C _i_ if + 1 / Sx o\ / Sx o\ 

9<P ~ °k BY 2 "A 2 1 ’ 2ct K Sy L x 0 J 4u 2 V 2 ) Uy ) 



(45) 



1 dx o Sx i 
2u 2 BY By 



2 2 
1 U e 



S 



x o K o 

2 



IV. Specification of u g (cp), ^(cp) 

To this point, whether a linearized or nonlinearized approach is 
taken, the solution can be developed analytically only providing u g , i 
are known functions of cp. We want to generate solutions for x(cpj'i'), K(cp,y)> 
hypothesize a dependence of upon K and unwind the transformation via the 
definition of cp in the following manner. 

dcp = v^. u g dx (46) 



J 1 v t (kX.Tj]) = I U e ^ 
which results implicitly in 



(47) 



f(cp) = g(x) 



After performing this integration, we can solve for cp(x) based on a given 
u g (x) and the hypothesized functional dependence for v (K) . The solutions 
can then be retransformed into physical variables and the validity of the 
hypotheses can be checked by comparison with the experimental u and K field. 
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In order to obtain expressions for and £ as functions of cp, we 
must first specify a general flow field to be studied. The following 
development serves only as a guide toward the specification of a class of 
functions which u e (c p) and j0 (cp) must be in order to adequately model 
realistic wake flows in pressure gradients. However, none of the approxi- 
mations in the following section involves an actual specification of v+ 
since u g and l could be specified as any arbitrary class of functions of 
cp and the actual physical flow field could later be inferred. 

The flow external to the wake in the experimental phase of the 
investigation was flowing through a channel whose area ratio is 



A_ = 1 

A 0 J 1 + px 



where 



P = 




( 48 ) 



( 49 ) 



For incompressible flow, this area distribution results in a velocity 
distribution 



— = J 1 + px (50) 

U 0 

with a resultant pressure gradient 

S - - *k> < 51) 

In order to map u g and i into the cp plane, we need a reliable approximate 
guide as to the shape of the v+ function so as to specify the transformation 

T> 
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from x to cp. From Schlichting (Ref 6) planar incompressible constant 
pressure wakes may be approximately solved with a of the form 

v t = k u g 9 (52) 

In this expression 9 serves as the characteristic length and for weak 
wakes is a measure of the macroscale {&) as is shown by the following 
development. From the definition of y we can write 

1 = / e I ( 1 ' fa) 2 dY (53) 

e 

which for weak wakes can be approximated by 

V “ I Te ( d + x) IT (54) 

0 x 2u ' 
e 

In addition from the definition of momentum thickness 



Y_ 



0 e e 0 e e 



(55) 



the following approximate expression can be derived 



u = e “C^ xdf 



(56) 



By comparing Equations 5^ and 56 



u A = Y q + u 9 
e e e 



(57) 



Since mass flux in the wake and hence Y g will be approximately conserved, 
9 will reflect the functional dependence of the macroscale on x. 
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In order to obtain an expression for 0 we can return to the momentum 

integral equation written with the specified external field which 

dp 

generates a constant pressure gradient ^ = - PIq. The momentum integral 
equation 



de 

dx 



n dU 

J e r 

U dx L H 



2 ] = 0 



(58) 



can be integrated for a constant shape factor (H) resulting in 



4 = (i ♦ M 1 ♦ h / 2 



(59) 



With this external velocity field, the Schlichting model for v+ ? and the 

"C 

solution for 0 we can write the transformation as 



x 



cp = k f u 0 dx 

J o e 



(6o) 



after inserting the expressions for 0 and u g we obtain 



k0 o u o 

9 = p(l - H/2) L (1 + 



0~0 f /, . - h /2 _ ■ 



(61) 



Using this relationship, the expressions for u g (x) and e(x) can be 
transformed into u g (cp) and 9 (cp) 



u e (cp) = u 0 (i - 






H 



— ) k Vo 



, 1/(2 

2 <P) 



- H) 



(62) 



e(cp) = e 0 ( x - 



(H + 2)/(H - 2) 



H 



— ) ke o u i 



2 V; 
0 



(63) 
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It should now be clear that we can consider the system of governing 
equations in the problem to be closed and well defined since all 
functional coefficients have been specified as functions of cp. Namely 
S(cp) and i(cp) 



SW = ^ = (1 ' (h H) ke o u o v) 



2/(H - 2) 



(64) 



and from Equation (57) 



Mcp)u e (cp) = Y e + u o e o (l - 



_p \ 


( 2 \ 


, 2 V 

k Vo 


U - 2. ) 



(H + 1)/(H - 2) 



(65) 



The specification of the value of the shape factor (H) to be used 
must be consistent with the formulation and the inherent approximations. 
From the definition of H we can write 




J ( x - $ * 



(66) 




(67) 



H = 




1 




e 



dy 




e 



(68) 
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Expanding the integrand, consistent with the order of the analysis we 
obtain 



H 




( 69 ) 



If we insert this value of H into Equations (6l), (64), and (65) we can 
write 



cp = 



2ke o u o 



[(1 + ex)2 - 1] 



(70) 



S(cp) (u ) i 1 + A P , 2 c p) 



2ke o u o 



(71) 



x(cp) u e (cp) = Y e + u Q e 0 (1 + 



2ke o u o 



cp) 



-2 



(72) 



SOLUTION OF THE EQUATIONS 

The general form for each of the linearized and nonli near i zed equa- 
tions for both K, x can te written 



H + f(cp)F + g(cp, y) 

9 dY 



(73) 



where f(cp) =0 for the momentum equations, 
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If we let F* = Fe"I f ^ cp ^ d£p , 
solution to the equation 



simple differentiation verifies that the 



ft . ill ♦ g {,, y)e -J f ^)*P ( 7 4) 

satisfies the general form of Equation (73) • 

The general solution of Equation (74) can he written in terms of 
Greens function G(x,y;cp,Y) as 



+co cp +oo 

F*(cp,y) = J* G(0,y;cp,y) I (y) dy + J J G(x,y;cp,Y) J (x,y) dx dy . (75) 



0 -o 



where the Greens function appropriate to the — - — — operator is 

Sep 3y 2 



G(x,y;cp,Y) = 1 exp I - .fo 7 v \ 

J 2rr(cp - x) L “ x) 



(76) 



Once the appropriate initial conditions (i) and nonhomogene ous functions (j) 
are inserted and integrated, the solution of F is 

F = Ffrel^ 9 ^ (77) 

COMPARISON WITH EXPERIMENT 

In order to test the adequacy with which the solutions presented in 
Equations (73) > (75 ) , and (77) predict a turbulent free mixing flow field, 
a sample quadrature was performed which could be compared with experiment. 
Figure 1 presents the initial conditions for x and K which were used in 
Equation (75) to solve for x 0 > Xj_> Kq, and - K 1 in cons ‘ tarrt adverse and 
favorable pressure gradients. Figures 1, 2, and 4 contain the results for 
a constant favorable pressure gradient with the experimental conditions in 
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indicated on the figures. Figures 1, 3> and 4 contain the results for 
the same calculations with a constant adverse pressure gradient of the 
same magnitude. The mathematical retransformation plots in both cases 
(Figures 4 and 5) were used with the eddy viscosity model 

v t ~ I ( 78 ) 

to test the ability of this model to reproduce experimental results for 
wakes in constant pressure gradients. The comparison of the analytical 
results with experiment shown in Figure 5 verifies that, with some further 
study of the proper values of the empirical constants, the approach out- 
lined here presents a consistent method for testing wake-type eddy 
viscosity models with streamwise pressure gradients and predicting 
untested physical situations with reliable eddy viscosity models. 

Although the analytical form of the solution has been obtained, 
often the initial conditions are discrete and do not satisfactorily fit 
any known analytic functions. In these cases a numerical solution of the 
governing equations can be easily obtained and the computer program 
needed to perform such a computation is listed and explained in the 
Appendix. 



CONCLUSIONS 

The equations of momentum and turbulent kinetic energy appropriate 
to free turbulent mixing have been developed. The resultant equations 
have been transformed into a plane which is independent of the eddy 
viscosity model and have been analytically solved by a perturbation 
technique. The solution depends upon prior knowledge of u e (cp) an( i J&(cp) 
and to specify these, an approximation for cp(x) must be available in 
order to study a flow field which is specified a priori. For this 
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analysis, the wake model given by Equation (52) proved to be satisfactory. 
If, on the other hand, only general classes of flow fields are of 
interest, u g (cp) and cp) may be specified independent of a specific 
physical problem with the resultant re-transformation determining the 
flow field which they implied. The results of the analysis have been 
compared to experimental wake data in constant adverse and favorable 
pressure gradients. The results of that comparison indicate that the 
analysis supplies a satisfactory method for obtaining analytical solu- 
tions to wake problems in freestream pressure gradients. 
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INITIAL CONDITIONS FOR SAMPLE CALCULATIONS. 
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SOLUTION FIELD FOR FAVORABLE PRESSURE 
GRADIENT. 
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FIG. 3 

SOLUTION FIELD FOR ADVERSE PRESSURE 
GRADIENT. - 
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1.5 -- 



1.0 — 



FAVORABLE PRESSURE GRADIENT 
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FIG. 4 

TRANSFORMATION PLOTS FROM <f> TO X. 
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(FEET) 



FIG. 5 

COMPARISON OF ANALYTICAL RESULTS 

WITH EXPERIMENT. 
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APPENDIX 

THE FCLLOWI NG COMPUTER PROGRAM IS A FORTRAN IV CODED FINITE DIFFERENCE 
CALCULATION OF THE ZEROTH AND FIRST ORDER SOLUTIONS TO THE EQUATIONS 
OF MOMENTUM AND TURBULENT KINETIC ENERGY APPLIED TO THE PROBLEM OF 
FREE TURBULENT MIXING IN FAVORABLE AND ADVERSE PRESSURE GRADIENTS 

TFF FOLLOWING PARAMETERS MUST BE INPUT AT THE BEGINNING OF THE DECK 
UINF U(1,J) EKl(l.J) SIGMA C Cl AL AO ALL MY MX DELY DELX D1 AKK 

:*c 7 +*#*** £ -$t ********■$ y* 



$ 5 $C 

* UINF IS THE INLET FREESTREAM VELOCITY * 

* * 

* UINF ALSO SERVES AS THE PERTURBATION UO * 

* $ 

* SIGMA IS THE CONSTANT IN ENERGY DIFFUSION * 

& ip. 

* C IS THE CONSTANT IN THE DISSIPATION TERM * 

* * 

* Cl IS THE CONSTANT IN THE EDDY VISCOSITY * 

* * 

* AL IS THE AREA OF THE TEST SECTION AT X=L * 

* AO IS THE AREA OF THE TEST SECTION AT X=0 * 

* ALL IS THE LENGTH OF THE TEST SECTION * 

ip 

* MY IS THE NUMBER OF LATERAL POINTS CALCULATED * 

afc ip 

* MX IS THE NUMBER OF AXIAL POINTS CALCULATED * 

* DELY IS THE LATERAL GRID SPACING * 

* D1 IS THE CONSTANT IN THE STEP SIZE * 

^ £ 

* DELX IS THE AXIAL SPACING TO RETRANSFORM * 

92 c ip 

* AKK IS THE CONSTANT IN THE APPROXIMATE PHI(X) * 

* * 

* MUMMY2 IS AN INDEX FOR A PARTICULAR PSI * 

afe * 



#***# 3 *#:*** a^***##*#*^#**##* ** 3* * 3k * * * ?lc * * # * * * >£ * 3fC # 



DIMENS ION 
DIMENSION 
DIMENSION 
DIMENSION 



CHI1 (200,81 ), CHI 2(200, 81), AK1( 200, 81), AK2( 200,81) 
THETA (200) ,DELSTR(200) ,S HAPE (200 ) , PS I E ( 200 ) , ELL ( 200 ) 
PHI ( 200) ,P2( 200) ,P4( 200) »X(200) ,FF(200) ,GG(2C0) 

APS I (81) ,AU(81),AK(81),Y(81),U(81),URATI0(81) ,EK1( 81) 
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DIMENSION PS I ( 8 1 ) 

C 

DO 1 1 = 1, MX 
DH 1 J = 1 « MY 

C 

CHI II I , J ) = 0. 0 
CH 1 2 ( I, J ) = 0.0 
AKKI , J) =0.0 
AK 2 ( I , J)=0.0 
C 

THETA ( I ) =0.0 
CEL STR ( I )= 0. 0 
S HAP E ( I )=0 .0 
PSIEU )=0.0 
ELL ( I )=0.0 
PHI ( I ) = 0 .0 
PSI ( I ) =0.0 
P? ( I ) = 0 . 0 
PA ( I )=0.0 
X( I ) = 0 . 0 
FF( I ) = 0.0 
GG( I )=0 .0 
C 

APSI ( J )=0.0 
AU ( J ) =0 . 0 
AK( J ) =0. 0 
Y ( J ) =0 . 0 
U ( J ) =0.0 
URATIO ( J )=0. 0 
EK1 (J)=0.0 
C 

1 CONTINUE 

C 

C THE FOLLOWING METHOD OF INPUTING U AND EK1 USES A ''STANDARD" SHAPE 
C FOR ROTH THE VELOCITY AND TURB ULE NT KINETIC ENERGY PROFILES ^OR USE 
C IN TEST CASES. IF THIS INPUT IS DESIRED ADDITIONAL INPUTS ARE NEEDED. 

C UCL-THE CENTERLINE VFLQCITY AT THE INLET, EKCL-THE TURBULENT KINETIC 
C FNFRGY CN THE CENTERLINE AT THE INLET, AND EKMAX-THE MAXIMUM TURBULENT 
C KINETIC ENERGY AT THE INLET. THESE VALUES SCALE THE SHAPE TO PROVIDE 
C SUITABLE INITIAL CONDITIONS FOR USE IN EVALUATING THE CONSTANTS AND 
C SPECIFYING A STEP SIZE WHICH WILL BE BOTH STABLE AND ACCURATE 
C 

U( 1) = 100.0 
U ( 2 ) =99 . 0 
U ( 3 ) =98. 0 
U l 4)=96. 0 
U (5 )=92 .0 
U( 6) =88.0 
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U ( 7) =81.0 
U( 8) = 73. 0 
U (9)=68 .0 
U( 10)=59.0 
U( 11 )=52.0 
U ( 12 ) =48 .0 
U( 13 ) =41 . 0 
U( 14 ) = 35 .0 
U { 1 5 ) =3 0 .0 
U( 16 ) =26. 0 
U( 17 ) = 2 1 .0 
U(18)=18.0 
U( 19)=14.0 
U ( 20 )= 1 1 .0 
U(21) =8.0 
U( 22) =5. 0 
U ( 23 ) =3 . 0 
U l 24) =2. 0 
U ( 25 ) = 1 . 0 
U (26 ) =0 . 5 
U ( 27) =0. 25 
U( 28)=0.0 

C 

ALPHA=1.0— (UCL/UINE) 

C 

nn 2 j = l f 5 1 

U( J) =UI.NF* (1.0-ALPHA*U (J)/100.0 ) 
2 CONTINUE 

C. 

FKl(l) =80.0 
EK 1 ( 2) = 8 1 . 0 
EK1 (3 ) =82 .0 
EK1 ( 4) =84. 0 
EK 1 ( 5 ) = 87. 0 
EK1 ( 6 ) =90 .0 
EK1( 7) =95. 0 
EK 1 ( 8 ) =98 . 0 
EK1 ( 9 ) =100 .0 
EK1 ( 10) =98.0 
EK1 ( 1 1 ) =92. 0 
EK1 (12) =83. 0 
EK1( 13 ) = 75. 0 
EK 1 ( 14 ) =67 . 0 
EK1( 15) =60.0 
EK1( 16)=53.0 
EK1 ( 1 7 ) =48 .0 
EK1 ( 18) =41.0 
EK 1 ( 1 9 ) =36 . 0 
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EK 1 ( 20 ) = 32 . 0 
EK1 121 )=26.0 
EK 1 ( 22) =21 .0 
EK 1 ( 23 ) = 17 . 0 
EK1 ( 24 ) = 13 .0 
EK1 ( 25) =10. 0 
EK 1 ( 26 ) = 8. 0 
EK1 (27 )=6.0 
EK1( 28) =4.0 
EK 1 ( 29 ) =2. 0 
EK1 (30 ) =1 .0 
EK1( 31 ) =0. 25 
EK1 ( 32 ) =0 *0 
C 

GAMMA=5.0*(1.0-(EKCL/EKMAX) ) 

00 344 J=l,9 

EK1 ( J) =EKMAX*' (1.0- (1. 0-EK1 ( J )/100 .0 )*GAMMA ) 

344 CONTINUE 

DO 345 J=10. 51 

EK1( J)=EKMAX*EK1 (JJ/10 0.0 

345 CONTINUE 
C 

UINF2=UINF**2 

C=C/C1 

BET A= ( (AO/AL )**2-l .0 )/ ALL 
C 

C SET UP INITIAL Y FIELD 

Y ( 1 ) =0 . 0 
DO 50 J =2 t MY 
Y(J)=FL0AT(J-1 )*DELY 

50 CONTINUE 
C 

C SET UP INITIAL PSI FIELD 

PSI (1 )=0.0 

PSI(2)=(U( 1)+U(2))*0.5*DELY 
DO 51 K = 3,MY 
K 1=K— 1 

A=0.5*( U( 1 ) +U( K ) ) 

B = 0 . 0 

DO 52 L =2 » K1 
B=U(L)+B 
52 CONTINUE 

PSI ( K ) = D E L Y‘‘ ( A + B ) 

51 CONTINUF 
C 

MY1=MY— 1 

DEL P SI = P S I ( MY) /MY1 
C 
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C SET UP UNIFORM PSI FIELD 

APSI (1 ) =0.0 
DO 53 J=2 ♦ MY 

APS I ( J ) =FLOAT( J-l ) *DELPSI 

53 CONTINUE 

r. 

WR IT E ( 6 . 99 5 ) 

995 FORMAT ( 1H1 t40X ♦ • «*** TABULATED INPUT VARIABLES IN PHYSICAL CCOP.DIN 
1ATES ***** ) 

C 

WRITc (6 .770 ) 

770 FORMAT! 1H0.25X, 'VELOCITY' ,10X, 'TURBULENT KINETIC ENERGY ', 10X ,' L AT E 
1RAL POSITION' . 10X, • STREAM FUNCTION') 

C 

00 880 J = 1 » MY 

WRITE (6, 996) U(Jl,EKl(J).Y(J) t APSI(J) 

996 FORMAT! IH0.25X , E15 . 7, 1 OX . E15 .7 , 10X , E15 . 7 ♦ 10X.E15.7) 

880 CONTINUF 

C 

C REDISTRIBUTE THE INPUT VALUES IN A UNIFORM PSI FIELD 

AU ( 1 ) =U ( 1 ) 

AK ( 1 ) = EK1 ( 1 ) 

DO 54 J = 2 » MY 
PC S= AP S I ( J ) 

AU ( J ) =P I F2 ( POS . P S UMY , U ) 

AK(J)=PIF2(P0S»PSI .MY.EK1) 

54 CONTINUE 
C 

C SET THE VALUES FOR THE DEPENDENT VARIABLES OF THE CALCULATION 

DO 55 J= 1 » MY 

CHIHl , J)=UINF**2-AU(J )**2 
AKK 1, J) =AK( J) 

PS I ( J ) = APS I ( J ) 

U ( J ) = AU ( J ) 

55 CONTINUE 
C 

DO 40 J = 1 « MY 
MUMMY= J 
MUMMY 1=J—1 

I F (CHI 1 (l.J) .LT.0.1 ) GO TO 41 

40 CONTINUF 

41 CONTINUE 
C 

PSIE ( 1 )=PSI (MUMMY) 

UE =U ( MUMMY ) 

UE2=UE**2 

C 

C SET UP NEW Y FIELD BASED ON UNIFROM PSI FIELD 
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Y ( 1 ) =0 . 0 

Y ( 2 J = ( ( U ( 1 ) +U ( 2 ) )/ (2.0*U(1 )*U{2 ) ) )*DELPS I 
DO 5 1 » K=3,MY 
Kl =K— 1 

A= 3. 5* ( ( 1 . 0/U( 1 ) )+(1.3/U(K))) 

B= 0 . 0 

DO 57 L = 2» Kl 
B=( 1.0/UtL) )+B 
57 CONTINUE 

Y(K)=DELPSI*(A+B) 

56 CONTINUE 
C 

ELL ( 1 ) =Y ( MUMMY ) 

C 

DO A3 J = 1 » M Y 

URATIOC J ) = SORT ( 1 . O-CHI 1( 1, J J/UE2) 

A3 CONTINUE 
C 

C CALCULATE DISPLACEMENT THICKNESS 

A = 0. 5*( ( 1. 3/URATI0( 1) )+(1.0/URATI0(MUMMY ) J-2.0) 

B=0 .0 

DC AA J=2 » MUMMY 1 
B = ( ( 1. 3/URATIOl J) )-l. 0 ) + B 
AA CONTINUE 

DFLSTR(1 )=DELPSI* (A+8J/UE 

r 

C CALCULATE MOMENTUM THICKNESS 

A = 0.5* (2.0-'JRATIC(l )-U RAT 10 ( MUMMY ) ) 

R = 0. 0 

DO A 5 J = 2» MUMMY 1 
B=(1.0— UPATIO(J) ) + B 
A 5 CONTINUE 

THETA( 1 ) = DELPS I* (A + BJ/UE 

0 

THETC=THETM 1) 

SHAPF(1)=DELSTR(1)/THETA(1 ) 

C 

WR IT E ( 6 . 99 7 ) 

997 FORMAT ( 1H1 » AOX » ' **** TABULATED INPUT VARIABLES IN STREAM FUNCTION 
1 COORD I NATE S ****' ) 

C 

WRITE (6,666) PH I ( 1 ) ,TH ET A ( I ), SHAPE ( I),DELSTR(1),ELL(1),PSIE(1) 

666 FORMAT( 1H0,2X» 'PHI = • , E 10. A , 2X , • THETA = • , E 10. A , 2X , • SHAP E FACTOR =• 
1 ,F10.A ,2X, • DELSTAR =' , E10.A,2X» 'MACROSCALE = • , E 1 0 . A, 2X , • MA SS FLUX 
2=' ,E10. A) 

C 

WP. ITE(6,771) 
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771 FORMATdHO, 20X , 'CHI', 15X, 'VELOCITY', 10X, 'TURBULENT KINETIC ENERGY' 
1 , 10X, • LATERAL POSITION' ,1 OX STREAM FUNCTION') 

C 

CO 881 J =1 » MY 

WRITE (6, 99 8) CHI1(1*J) ♦ AU ( J ) , AK1 ( 1 , J ) , Y ( J ) , PS I ( J ) 

998 FORMAT! 1HO, 10X , E 1 5 . 7 , 1 OX , E 1 5 . 7 , 10X , E 15. 7 , 1 OX , El 5. 7 , 1 OX , E 1 5 . 7 ) 

881 CONTINUE 
C 

C SET STABLE STEP SIZE FOR DELPHI 

DELPHI =D1* ( DEL PS 1**2) 

C 

PHI ( 1 ) =0 .0 
CO 5 1=2, MX 

PHI ( I ) =F LOAT(I -1)*DELPHI 
5 CONTINUE 
C 

AMEGA=0. 5*BETA/( AKK*THETO* (UINF**2) ) 

C 

CO 7 1=1 , MX 

P2(I ) = 1.0/( ( 1, 0+ A MEGA* PH 1(1) )**2) 

P4( I )= 1 .0/ ( ( 1 • 0+ AME GA* PH I ( I ) )** 4 ) 

7 CONTINUE 

0 

MX 1=MX— 1 
C 

C START MAIN CALCULATION LOOP 

DO 20 1=1, MX1 
C 

F =— C / ( ( P S I E ( I)+UINF*THET0*P2 (I ) )**2) 

C 

DO 10 J=2 , MY1 

DER I V= ( CHI 1 ( I , J + 1) + CHI1(I , J-l ) -2.0*CHI 1 ( I, J) )/ ( DEL PS 1**2 ) 

CHI 1 (1+1, J )=CHI1( I,J)+DELPHI*DERIV 

10 CONTINUE 
C 

DERIV=2.0*(CHI 1 ( 1,2) — C HI 1( I» 1) ) / ( DELPSI ** 2 ) 

CHI1 ( 1+1 »1 )=CHI1 ( I , 1)+DELPHI*DERIV 

C 

DO 11 J=2, MY1 

G=— 0 . 5*P2 ( I )*CHI 1( I ,J )*(CHI 1( I« J+1)-2.0*CHI1( I , J )+CHIl( I , J-l) )/ 
1(DELPSI**2) 

DER IV= ( C HI 2 ( I , J +1 )-2. 0*CHI 2( I , J ) +CH 1 2 ( I , J-l ) ) /( DE LPSI** 2 ) 

CHI2H+1 ,J)=CHI2 ( I ,J )+DELPHI* (DERIV+G) 

11 CONTINUE 
C 

G=— 0. 5 * ■ P2 ( I )*CH 1 1 ( I ,1 )*2.0*(CHI1 { I*2)-CHI 1( 1 , 1 ) )/ ( DELPS I **2 ) 
DERIV=2.0*(CHI 2 ( I , 2 )-CHI 2( I ♦ 1) ) / ( DELPSI **2 ) 

CHI2 (1+1,1 )=CH 12 ( I , 1 1+DELPH I *( DERIV+G) 
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c 



c 

c 



c 



c 

c 



c 



c 



DO 12 J =2 , MY1 

G=0.25*P2( I)*( ( (CHI 1(1 ,J+1)-CHI1 (I , J-l) ) / ( 2.0*DELPSI*UINF) )**2) 
CERIV=( AK1 ( I,J+1 )-2.0*AKl( I, J)+AK1( I * J— 1 ) ) /(DELPSI**2) 

AK1(I + 1,J) =AK1 ( I , J )+DELPHI*( ( DER I V/ S I GM A ) + F* AK1 ( I * J )+G) 

12 CONTINUE 

DFRIV=2.0*( AK1 (1*2) -AK1 (1,1) )/ ( DELPS 1**2 ) 

AKKI+l, 1) =AK1( 1,1 ) +DE LPH I * ( ( DER I V/ S I GMA ) +F* AK1 ( I , 1 ) ) 

DO 13 J=2 , MY1 

G1A=CHI1(I,J)=MAK1( I,J + 1)-2.0*AK1(I,J) + AK1(I ,J-1) )/(DELPSI**2) 

G1 B=( CHI1 ( I, J + l ) — CH III 1 , J-l) )*( AK1( I , J+1J-AKH l, J-l) )/f 4.0* 

1 (DELPSI**2) ) 

Gl=-0. 5*P2 ( I )* ( G 1A+G1B ) 

G2 =0 . 1 2 5*P4 ( I)*CHI1(I,J)*(((CHI1(I,J+1I-CHI1( I , J- 1 ) ) /( 2 . 0*OE LP SI* 

1 U INF ) )**2) 

G3=0.5*P2( I )*(CHI 1( I » J + 1 ) —CHI 1 ( I , J-l ) )* ( CHI 2 ( I , J+l) -CHI 2 ( I , J-l ) )/ 

1 (4.0* (DELPS 1**2 )*(U IN F**2 ) ) 

G4=- 0 . 5*C*P4( I )*CHI 1(1 » J ) * AK1 ( I,J)/((PSIE ( I ) +UI NF*THET0*P2 ( I ) )**2 ) 
G=G1+G2+G3+G4 

DERIV=( AK2 ( I , J +1 )-2 .3* AK2 ( I , JJ+AK2 ( I , J— 1 ) )/( DELPS 1**2 ) 

AK2 ( 1 + 1, J)=AK2( I » J J+DELPHI* ( ( DE R I V/S I GMA ) + F* AK2 ( I , J ) +G ) 

13 CONTINUE 

G 1A=CH I 1(I,1)*2.0*(AK1(I,2)-AK1(I,1) ) / ( DE L PS 1**2 ) 



G1 FS=0 . 0 

Gl=-0. 5*P2( I )* ( G1A+G1B) 

G2=0 . 0 
G3 =0 . 0 

G4=-0. 5#C* P4 ( I )*CHI 1(1 ,1)*AK1( I ,1 )/( (PSIE( I )+U INF+THET0+P2 ( I ) )**2) 
G=G1+G2+G3+G4 

CER I V=2 .0* ( AK2 (1,2 )-AK2( 1,1) )/ ( DELPS 1**2 ) 

AK2( 1+1,1) =AK2 ( I ,1)+DELPHI* ( ( DE R I V/S I GMA ) * F* AK2 ( 1,1 )+G) 

END OF MAIN CALCULATION LOOP 



DO 75 J = 1 « MY 

U( J) = S0RT(UINF2/P2( I ) — CHI 1 ( 1 + 1 , J ) ) 

AU(J)=SQRT (UINF2/P2 (I)-(CHIl ( 1+1 , J ) +CHI2 ( I +1 , J )/U INF2 ) ) 
75 CONTINUE 



DO 61 J=1,MY 
MUMMY= J 
MUMMY 1 = J— 1 

IF(CHIl(I+l,J).LT.D.l) GO TO 70 
61 CONTINUE 



70 CONTINUE 
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MUMMY2=MUMMY/2 

PSIE (1+1 )=PSI ( MUMMY ) 

. U C =U ( MUMMY ) 

UE2=UE**2 

CALCULATE NEW Y FIELD 
Y ( 1 ) =0 . 0 

Y ( 2) =0. 5*DELPSI* ( 1 . 0/U ( 1 ) + 1 . 0/U ( 2 ) ) 

DO 59 K = 3, MY 
K 1=K-1 

A=0.5*< 1.0/U( 1) + 1. 0/U( K) ) 

B=0 . 0 

DO 60 L =2 « K1 
B= 1 • 0/U ( L ) +B 
60 CONTINUE 

Y(K) =DE L PS I* (A + B) 

59 CONTINUE 

ELL ( I + 1 ) =Y ( MUM MY ) 

DO 62 J=l, MY 

URATIO( J)=SORT <1 .0-CHI1 (1+1, J )/UE2) 

APSI(J)=SQRT(1. 0— (CHI 1 ( 1 + 1 , J )+CHI 2 ( I + 1 , J ) /UI NF2 ) / UE2 ) 

EK1 ( J ) = AK1 ( I +1 , J J+AK2I I + l f J J/UINF2 

62 CONTINUE 

A=0 .5* ( ( 1.0/UR AT 10(1) ) + ( 1.0/UR AT 10 (MUMMY) )-2. 0) 

B = 0. 0 

DO 63 J = 2, MUMMY1 
B = ( ( 1 • O/URAT 10 ( J ) ) - 1 .0 )+B 

63 CONTINUE 

DELSTR( I+1)=DELPSI*(A+B)/UE 

A=0 . 5* ( 2 . 0— URATI 0( 1)-URATI0( MUMMY ) ) 

B=0 .0 

DO 64 J = 2 » MUMMY1 
B= ( 1. O-URATIO ( J ) ) + B 

64 CONTINUE 

THETA< 1+1) =DEL PS I* (A+BI/UE 

SHAPE ( I+1)=DELSTR( 1+1) /THETA ( I +1) 

C 

C 

WRITE( 6,991 ) 

991 FORMAT ( 1H1, 40X, ***** INTEGRATED PROFILE PARAMETERS **** * ) 
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WRITE (6, 992) P H I ( I +1 ) , T HET A( I + 1 ) , SHAPE! I + 1 ) , DEL STR ( I + 1 ) , ELL ( I + 1 ) , 
IPS 15(1 + 1) 

992 FORMAT ( 1 HO ,2X, ' PHI = • , E10.4, 2X, 'THETA = • , E 10. 4, 2X , • SHAP E FACTOR =• 
1 ,E10.4,2X,' DELSTAR =• , E 10. 4 , 2X , ' MACROS CAL E = • , El 0 .4, 2X, • M ASS FLUX 
2=' , E10.4) 

C 

WRITE( 6,999) 

999 FORMAT ( 1H0, 40X , • if*** CALCULATED OUTPUT VARIABLES **+*• 1 
C 

WRITE (6,772) 

772 FORMAT ( 1H0, » J' , 8X, 'U/UE • ,8X, ' UT/UE' ,8X,' U' , 8 X , ' UT ' , 8X , • Y • ,8X , • PS I ' 
1 ,8X, • CHI1' ,8X, 'CHI2"«8X« 'AK1N8X, «AK2', 8X, 'AKT* ) 

C 

DO 883 J=1,MY 

WRITE (6, 990) J»URATIO(J),APSI(J),U(J),AU(J),Y(J),PSI(J) ,CHI 1( I+1,J 
1),CHI2(I+1,J),AK1(I+1,J) , AK2 ( I + 1 » J ) , EK1 ( J ) 

9 90 FORMAT ( 1 HO , I2,2X,F5.3,2X,F5.3,2X,F7.3,2X,F7.3,2X,F5.3,2X,F7.3,2X,E 
112.5,2X»E12.5,2X,E12.5,2X,E12.5,2X»E12.5) 

883 CONTINUE 
I C ALC= I 

I F ( U ( 1 ) .LT.5.0) GO TO 21 

C 

20 CONTINUE 
C 

21 CONTINUE 
WRITE (6,993) 

993 FORMAT( 1H1 , ***** SUMMARY OF INTEGRATED QUANTITIES ****') 

C 

DO 874 I =1 , I CA LC 

WR I TE ( 6,994) THE TA ( I ) , DELSTR( I ) .SHAPE (I ) , PS I E ( I ) , ELL ( I ) , I 

994 FORMAT ( 1H0, 5E20.8, 13) 

874 CONTINUE 

C 

X( 1 ) = 0.0 

DC 869 I =2, MX 

X( I )=FLOAT( I-D + DELX 

869 CONTINUE 
C 

FF( 1 ) =0. 0 
DO 870 1=2, MX 

FF ( I )=2.0*C1*U INF* ( (1 .0+BETA*X ( I ) )**1 .5- 1 . 0 )/ ( 3 . 0*B ETA) 

870 CONTINUE 
C. 

GG ( 1 )=0.0 

GG ( 2 )= 0. 5* ( 1. 0/(S0RT(AKl( 1,MUMMY2) )*ELL(1) )+ 1.0/ (SORT ( AK1 (2 ,MUMM 

1Y2))*ELL(2)) )*DELPHI 
DO 871 I =3 , ICA LC 
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1 1=1-1 

A = 0. 5*( 1.0/ ( SQRKAKK 1 ,MUMMY2 ) ) *E LL ( 1 ) ) + 1.0/ (SORT ( AK1 ( I , MUM MY 2 ) )* 
1 ELL C I ) ) ) 

9 = 0.0 

DO 872 K = 1 , I 1 

8=1.0/ (SORT (AK1 ( I,MUMMY2))*ELL(K) ) + B 

872 CONTINUE 

GG ( I )=DELPHI*(A+B) 

871 CONTINUE 

DO 877 I =1 C ALC , MX 
GG( I )=0.0 
877 CONTINUE 

WRITE (6. 850 ) Cl 

850 FORMAT ( 1H1 , 30X , ' RETRANSFORMATION RESULTS FOR Cl =',E10.3) 

W R IT E ( 6 . 848) 

848 FORMAT (1 HO, 30X, «F» , 20X,'G' , 20X, 'PHI' , 20X, 'X' ) 

DO 849 1=1, MX 

WRITE (6 ,873 ) FF(I),GG(I),PHI(I),X(I) 

873 FORMAT( 1H0, 10X ,4E20.8) 

849 CONTINUE 

END 



FUNCTION PIF2 ( X , X L 1ST , N , FL I ST ) 

C FUNCTION P I F 2 IS A SECOND ORDER LOOKUP FUNCTION 
DIMENSION XL 1ST (ICO), FLIST (100) 

BLIF (P,Q,R,S,T) = ( (Q-P)* ( S— T ) / (R-Q)+S ) 

IF ( X-XLIST(N) ) 2,1,1 

1 I = Nt- 1 
GO TO 5 

2 I F ( X— XL I ST ( 1 ) ) 4,4,6 

4 1=1 

5 K = 1 

GO TO 30 

6 K = 2 

7 DO 8 I = 1,N 

IF (X-XLISTI I ) ) 9,9,8 

8 CONTINUE 
I =N 

91=1-1 

30 BLIF1 = BLIF(X,XLIST(I ),XLIST( I+1),FLIST( I ) , FL 1ST ( 1+ 1 ) ) 

10 IF ( K— 1 ) 11,11,12 

11 PIF2 = BLIF1 
RETURN 
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12 I F ( ( I +2 ) -N ) 13,13,16 

13 IF (II-l)-l) 15,14,14 

14 I F ( ABS ( XLIST ( I — 1 )-X)-ABS(XLI ST ( I+2J-X) ) 16, 15,15 

15 L = 1+2 
GO TO 17 

16- L = 1-1 

17 BLIF2 = BLIF ( X , XL I ST ( I ) ,XL 1ST ( L ) , FL 1ST ( I ) , FL 1ST ( L ) ) 
PIF2 = BLIF ( X, XLI ST( 1+1) ,XLIST(L) , B LI F 1 , B L I F2 ) 

18 RETURN 
END 
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